Question 1

Use R package UScensus2010county to complete the following tasks: (20 pt.)

Question 1(a)

Plot a map of New York counties using the plot function.

if (!require(UScensus2010)) {install.packages("UScensus2010")}
library(UScensus2010)
if (!require(UScensus2010county)) {install.county("osx")}
if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010county)
library(UScensus2010tract)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df 
plot(shp)

Question 1(b)

Plot a map of New York counties using the qtm function.

library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df
plot(shp)

if (!require(tmap)) {install.packages("tmap")}
library(tmap)
qtm(shp)

Question 1(c)

How many counties in New York State?

shp <- new_york.county10
df <- shp@data
nrow(df)
[1] 62

Question 1(d)

What’s the 3-digit fips code of Broome County?

shp <- new_york.county10
df <- shp@data
sel <- df$fips[df$NAME10 == "Broome"]
sel
[1] "36007"
sprintf('%03d', as.numeric(sel) %% 1000)
[1] "007"

Question 1(e)

Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.

library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
pop <- df$P0010001
sum(pop)
[1] 19378102
min(pop)
[1] 4836
max(pop)
[1] 2504700
mean(pop)
[1] 312550
median(pop)
[1] 91301
library(e1071)
skewness(pop)
[1] 2.517639

Question 1(f)

Create a histogram and a boxplot of the population.

hist(pop)

boxplot(pop)

Question 2

Use R package UScensus2010tract to complete the following tasks: (20 pt.)

if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010tract)

Question 2(a)

Plot a map of New York census tracts using the plot function.

data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
df 
plot(shp)

Question 2(b)

Compute the total population based on census tracts.

data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
pop <- sum(df$P0010001)
pop
[1] 19378102

Question 2(c)

Select all census tracts in Broome County and plot the map.

data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
sel <- df$county == "007"
plot(shp[sel,])

Question 2(d)

What’s the total population of Broome County?

Broomepop <- df$P0010001[df$county == "007"]
sum(Broomepop)
[1] 200600

Question 2(e)

Create a histogram and a boxplot of population based on census tracts of Broome County.

hist(Broomepop)

boxplot(Broomepop)

Question 2(f)

Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.

dfnew <- df[,1:5]
shp@data <- dfnew
shp@data
shp$ratio <- shp$P0010001/sum(Broomepop)
writeOGR(shp,dsn = "output",layer = "newshp",driver = "ESRI Shapefile",overwrite_layer = TRUE)
NAs introduced by coercion
newshp <- readOGR(dsn = "output", layer = "newshp")
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output", layer: "newshp"
with 4919 features
It has 6 fields
newshp@data

Question 3

Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)

if (!require(ncdf4)) {install.packages("ncdf4")}
library(ncdf4)
library(raster)

Question 3(a)

Load the multi-band raster dataset “NDVI.nc” into R.

ndvi <- brick("NDVI.nc")
ndvi
class       : RasterBrick 
dimensions  : 360, 720, 259200, 36  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\NDVI.nc 
names       : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ... 
Date        : 2000-01-01, 2002-12-01 (min, max)
varname     : NDVI 

Question 3(b)

Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.

nrow(ndvi)
[1] 360
ncol(ndvi)
[1] 720
ncell(ndvi)
[1] 259200
extent(ndvi)
class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -90 
ymax        : 90 
bbox(ndvi)
    min max
s1 -180 180
s2  -90  90
res(ndvi)
[1] 0.5 0.5
projection(ndvi)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Question 3(c)

Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.

ndvimean <- mean(ndvi)
ndvimean
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.8537528  (min, max)
ndvimax <- max(ndvi)
ndvimax
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.9681  (min, max)
writeRaster(ndvimean, "output/ndvimean.nc", overwrite = TRUE) 
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output\ndvimean.nc 
names       : layer 
zvar        : layer 
writeRaster(ndvimax, "output/ndvimax.nc", overwrite = TRUE) 
class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output\ndvimax.nc 
names       : layer 
zvar        : layer 

Question 3(d)

Plot the maps of monthly NDVI of the year 2001

ndvi2001 <- ndvi[[13:24]] 
plot(ndvi2001)

Question 3(e)

Create histograms of monthly NDVI of the year 2001.

hist(ndvi2001)
size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

Question 3(f)

Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.

plot(ndvi, 7)
values <- click(ndvi, n=1, xy=FALSE)

plot(as.vector(values), type = "b")
no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfError in plot.window(...) : need finite 'xlim' values

LS0tDQp0aXRsZTogIkdFT0ctNTMzIExhYiA5Ig0KYXV0aG9yOiBZdXJ1aSBaaGFuZw0Kb3V0cHV0OiANCiAgaHRtbF9ub3RlYm9vazoNCiAgICB0b2M6IFRSVUUNCiAgICB0b2NfZmxvYXQ6IFRSVUUNCi0tLQ0KDQojIyBRdWVzdGlvbiAxDQpVc2UgUiBwYWNrYWdlICoqVVNjZW5zdXMyMDEwY291bnR5KiogdG8gY29tcGxldGUgdGhlIGZvbGxvd2luZyB0YXNrczogICgyMCBwdC4pDQoNCiMjIyBRdWVzdGlvbiAxKGEpDQpQbG90IGEgbWFwIG9mIE5ldyBZb3JrIGNvdW50aWVzIHVzaW5nIHRoZSAqKnBsb3QqKiBmdW5jdGlvbi4NCmBgYHtyfQ0KaWYgKCFyZXF1aXJlKFVTY2Vuc3VzMjAxMCkpIHtpbnN0YWxsLnBhY2thZ2VzKCJVU2NlbnN1czIwMTAiKX0NCmxpYnJhcnkoVVNjZW5zdXMyMDEwKQ0KaWYgKCFyZXF1aXJlKFVTY2Vuc3VzMjAxMGNvdW50eSkpIHtpbnN0YWxsLmNvdW50eSgib3N4Iil9DQppZiAoIXJlcXVpcmUoVVNjZW5zdXMyMDEwdHJhY3QpKSB7aW5zdGFsbC50cmFjdCgib3N4Iil9DQpsaWJyYXJ5KFVTY2Vuc3VzMjAxMGNvdW50eSkNCmxpYnJhcnkoVVNjZW5zdXMyMDEwdHJhY3QpDQpkYXRhKG5ld195b3JrLmNvdW50eTEwKQ0Kc2hwIDwtIG5ld195b3JrLmNvdW50eTEwDQpkZiA8LSBzaHBAZGF0YQ0KZGYgDQpwbG90KHNocCkNCg0KYGBgDQoNCiMjIyBRdWVzdGlvbiAxKGIpCQ0KUGxvdCBhIG1hcCBvZiBOZXcgWW9yayBjb3VudGllcyB1c2luZyB0aGUgKipxdG0qKiBmdW5jdGlvbi4NCmBgYHtyfQ0KbGlicmFyeShVU2NlbnN1czIwMTBjb3VudHkpDQpkYXRhKG5ld195b3JrLmNvdW50eTEwKQ0Kc2hwIDwtIG5ld195b3JrLmNvdW50eTEwDQpkZiA8LSBzaHBAZGF0YQ0KZGYNCnBsb3Qoc2hwKQ0KaWYgKCFyZXF1aXJlKHRtYXApKSB7aW5zdGFsbC5wYWNrYWdlcygidG1hcCIpfQ0KbGlicmFyeSh0bWFwKQ0KcXRtKHNocCkNCg0KYGBgDQoNCg0KIyMjIFF1ZXN0aW9uIDEoYykJDQpIb3cgbWFueSBjb3VudGllcyBpbiBOZXcgWW9yayBTdGF0ZT8NCmBgYHtyfQ0Kc2hwIDwtIG5ld195b3JrLmNvdW50eTEwDQpkZiA8LSBzaHBAZGF0YQ0KbnJvdyhkZikNCmBgYA0KDQojIyMgUXVlc3Rpb24gMShkKQkNCldoYXTigJlzIHRoZSAzLWRpZ2l0ICoqZmlwcyoqIGNvZGUgb2YgQnJvb21lIENvdW50eT8NCmBgYHtyfQ0Kc2hwIDwtIG5ld195b3JrLmNvdW50eTEwDQpkZiA8LSBzaHBAZGF0YQ0Kc2VsIDwtIGRmJGZpcHNbZGYkTkFNRTEwID09ICJCcm9vbWUiXQ0Kc2VsDQoNCnNwcmludGYoJyUwM2QnLCBhcy5udW1lcmljKHNlbCkgJSUgMTAwMCkNCg0KYGBgDQoNCiMjIyBRdWVzdGlvbiAxKGUpCQ0KQ29tcHV0ZSBkZXNjcmlwdGl2ZSBzdGF0aXN0aWNzIG9mIHRoZSBwb3B1bGF0aW9uIGNvbHVtbiAoKipQMDAxMDAwMSoqKSwgaW5jbHVkaW5nIHRvdGFsLCBtaW5pbXVtLCBtYXhpbXVtLCBtZWFuLCBtZWRpYW4sIGFuZCBza2V3bmVzcy4gDQpgYGB7cn0NCmxpYnJhcnkoVVNjZW5zdXMyMDEwY291bnR5KQ0KZGF0YShuZXdfeW9yay5jb3VudHkxMCkNCnNocCA8LSBuZXdfeW9yay5jb3VudHkxMA0KZGYgPC0gc2hwQGRhdGENCnBvcCA8LSBkZiRQMDAxMDAwMQ0Kc3VtKHBvcCkNCm1pbihwb3ApDQptYXgocG9wKQ0KbWVhbihwb3ApDQptZWRpYW4ocG9wKQ0KDQpsaWJyYXJ5KGUxMDcxKQ0Kc2tld25lc3MocG9wKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAxKGYpCQ0KQ3JlYXRlIGEgaGlzdG9ncmFtIGFuZCBhIGJveHBsb3Qgb2YgdGhlIHBvcHVsYXRpb24uDQpgYGB7cn0NCmhpc3QocG9wKQ0KYm94cGxvdChwb3ApDQpgYGANCg0KDQojIyBRdWVzdGlvbiAyDQpVc2UgUiBwYWNrYWdlICoqVVNjZW5zdXMyMDEwdHJhY3QqKiB0byBjb21wbGV0ZSB0aGUgZm9sbG93aW5nIHRhc2tzOiAgICAoMjAgcHQuKQ0KYGBge3J9DQppZiAoIXJlcXVpcmUoVVNjZW5zdXMyMDEwdHJhY3QpKSB7aW5zdGFsbC50cmFjdCgib3N4Iil9DQpsaWJyYXJ5KFVTY2Vuc3VzMjAxMHRyYWN0KQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAyKGEpCQ0KUGxvdCBhIG1hcCBvZiBOZXcgWW9yayBjZW5zdXMgdHJhY3RzIHVzaW5nIHRoZSAqKnBsb3QqKiBmdW5jdGlvbi4NCmBgYHtyfQ0KZGF0YSgibmV3X3lvcmsudHJhY3QxMCIpDQpzaHAgPC0gbmV3X3lvcmsudHJhY3QxMA0KZGYgPC0gc2hwQGRhdGENCmRmIA0KcGxvdChzaHApDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDIoYikNCkNvbXB1dGUgdGhlIHRvdGFsIHBvcHVsYXRpb24gYmFzZWQgb24gY2Vuc3VzIHRyYWN0cy4NCmBgYHtyfQ0KZGF0YSgibmV3X3lvcmsudHJhY3QxMCIpDQpzaHAgPC0gbmV3X3lvcmsudHJhY3QxMA0KZGYgPC0gc2hwQGRhdGENCnBvcCA8LSBzdW0oZGYkUDAwMTAwMDEpDQpwb3ANCmBgYA0KDQojIyMgUXVlc3Rpb24gMihjKQ0KU2VsZWN0IGFsbCBjZW5zdXMgdHJhY3RzIGluIEJyb29tZSBDb3VudHkgYW5kIHBsb3QgdGhlIG1hcC4gDQpgYGB7cn0NCmRhdGEoIm5ld195b3JrLnRyYWN0MTAiKQ0Kc2hwIDwtIG5ld195b3JrLnRyYWN0MTANCmRmIDwtIHNocEBkYXRhDQpzZWwgPC0gZGYkY291bnR5ID09ICIwMDciDQpwbG90KHNocFtzZWwsXSkNCmBgYA0KDQojIyMgUXVlc3Rpb24gMihkKQ0KV2hhdOKAmXMgdGhlIHRvdGFsIHBvcHVsYXRpb24gb2YgQnJvb21lIENvdW50eT8NCmBgYHtyfQ0KQnJvb21lcG9wIDwtIGRmJFAwMDEwMDAxW2RmJGNvdW50eSA9PSAiMDA3Il0NCnN1bShCcm9vbWVwb3ApDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDIoZSkNCkNyZWF0ZSBhIGhpc3RvZ3JhbSBhbmQgYSBib3hwbG90IG9mIHBvcHVsYXRpb24gYmFzZWQgb24gY2Vuc3VzIHRyYWN0cyBvZiBCcm9vbWUgQ291bnR5Lg0KYGBge3J9DQpoaXN0KEJyb29tZXBvcCkNCmJveHBsb3QoQnJvb21lcG9wKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAyKGYpDQpTZWxlY3QgdGhlIGZpcnN0IGZpdmUgY29sdW1ucyBvZiB0aGUgc2hhcGVmaWxlIG9mIEJyb29tZSBDb3VudHkgY2Vuc3VzIHRyYWN0OyBhZGQgYSBuZXcgcG9wdWxhdGlvbiByYXRpbyBjb2x1bW4gKD0gY2Vuc3VzIHRyYWN0IHBvcHVsYXRpb24gLyBjb3VudHkgcG9wdWxhdGlvbik7IHNhdmUgdGhlIG5ldyBzaGFwZWZpbGUgdG8gdGhlIGhhcmQgZHJpdmUuIA0KYGBge3J9DQpkZm5ldyA8LSBkZlssMTo1XQ0Kc2hwQGRhdGEgPC0gZGZuZXcNCnNocEBkYXRhDQpzaHAkcmF0aW8gPC0gc2hwJFAwMDEwMDAxL3N1bShCcm9vbWVwb3ApDQoNCndyaXRlT0dSKHNocCxkc24gPSAib3V0cHV0IixsYXllciA9ICJuZXdzaHAiLGRyaXZlciA9ICJFU1JJIFNoYXBlZmlsZSIsb3ZlcndyaXRlX2xheWVyID0gVFJVRSkNCm5ld3NocCA8LSByZWFkT0dSKGRzbiA9ICJvdXRwdXQiLCBsYXllciA9ICJuZXdzaHAiKQ0KbmV3c2hwQGRhdGENCmBgYA0KDQoNCiMjIFF1ZXN0aW9uIDMNCg0KVXNlIFIgcGFja2FnZXMgKipyYXN0ZXIqKiBhbmQgKipuY2RmNCoqIHRvIGNvbXBsZXRlIHRoZSBmb2xsb3dpbmcgdGFza3M6ICAgICAoMjAgcHQuKQ0KYGBge3J9DQppZiAoIXJlcXVpcmUobmNkZjQpKSB7aW5zdGFsbC5wYWNrYWdlcygibmNkZjQiKX0NCmxpYnJhcnkobmNkZjQpDQpsaWJyYXJ5KHJhc3RlcikNCmBgYA0KDQojIyMgUXVlc3Rpb24gMyhhKQkJDQpMb2FkIHRoZSBtdWx0aS1iYW5kIHJhc3RlciBkYXRhc2V0IOKAnE5EVkkubmPigJ0gaW50byBSLg0KYGBge3J9DQpuZHZpIDwtIGJyaWNrKCJORFZJLm5jIikNCm5kdmkNCmBgYA0KDQojIyMgUXVlc3Rpb24gMyhiKQkJDQpHZXQgdGhlIGJhc2ljIGluZm9ybWF0aW9uIGFib3V0IHRoZSBkYXRhc2V0LCBpbmNsdWRpbmcgdGhlIG51bWJlciBvZiByb3dzLCBjb2x1bW5zLCBjZWxscywgYW5kIGJhbmRzOyBzcGF0aWFsIHJlc29sdXRpb24sIGV4dGVudCwgYm91bmRpbmcgYm94LCBhbmQgcHJvamVjdGlvbi4NCmBgYHtyfQ0KbnJvdyhuZHZpKQ0KbmNvbChuZHZpKQ0KbmNlbGwobmR2aSkNCmV4dGVudChuZHZpKQ0KYmJveChuZHZpKQ0KcmVzKG5kdmkpDQpwcm9qZWN0aW9uKG5kdmkpDQpgYGANCg0KIyMjIFF1ZXN0aW9uIDMoYykJDQpBZ2dyZWdhdGUgYWxsIGJhbmRzIHRvIGdlbmVyYXRlIGEgbWVhbiBORFZJIHJhc3RlciBhbmQgYSBtYXhpbXVtIE5EVkkgcmFzdGVyOyBzYXZlIHRoZSB0d28gbmV3IHJhc3RlciBkYXRhc2V0cyB0byB0aGUgaGFyZCBkcml2ZS4gDQpgYGB7cn0NCm5kdmltZWFuIDwtIG1lYW4obmR2aSkNCm5kdmltZWFuDQoNCm5kdmltYXggPC0gbWF4KG5kdmkpDQpuZHZpbWF4DQoNCndyaXRlUmFzdGVyKG5kdmltZWFuLCAib3V0cHV0L25kdmltZWFuLm5jIiwgb3ZlcndyaXRlID0gVFJVRSkgDQp3cml0ZVJhc3RlcihuZHZpbWF4LCAib3V0cHV0L25kdmltYXgubmMiLCBvdmVyd3JpdGUgPSBUUlVFKSANCmBgYA0KDQojIyMgUXVlc3Rpb24gMyhkKQkJCQ0KUGxvdCB0aGUgbWFwcyBvZiBtb250aGx5IE5EVkkgb2YgdGhlIHllYXIgMjAwMQ0KYGBge3J9DQpuZHZpMjAwMSA8LSBuZHZpW1sxMzoyNF1dIA0KcGxvdChuZHZpMjAwMSkNCmBgYA0KDQojIyMgUXVlc3Rpb24gMyhlKQkNCkNyZWF0ZSBoaXN0b2dyYW1zIG9mIG1vbnRobHkgTkRWSSBvZiB0aGUgeWVhciAyMDAxLg0KYGBge3J9DQpoaXN0KG5kdmkyMDAxKQ0KYGBgDQoNCiMjIyBRdWVzdGlvbiAzKGYpCQkJDQpQbG90IHRoZSBORFZJIG1hcCBvZiBKdWx5IDIwMDA7IGNsaWNrIGFueSBsb2NhdGlvbiB3aXRoIGRhdGEgb24gdGhlIG1hcCBhbmQgcmV0cmlldmUgdGhlIE5EVkkgdGltZSBzZXJpZXMgZm9yIGFsbCB5ZWFyczsgcGxvdCB0aGUgTkRWSSB0aW1lIHNlcmllcyBvZiB0aGUgc2VsZWN0ZWQgbG9jYXRpb24uIA0KYGBge3J9DQpwbG90KG5kdmksIDcpDQp2YWx1ZXMgPC0gY2xpY2sobmR2aSwgbj0xLCB4eT1GQUxTRSkNCnBsb3QoYXMudmVjdG9yKHZhbHVlcyksIHR5cGUgPSAiYiIpDQpgYGANCg0K